### Load packages
library(tidyverse)
library(readxl)
require(forecast)
library(zoo)
library(kableExtra)
library(astsa)
### Load data
gved_y_raw <- read_excel("data/nipa_3155.xlsx")
gv_y_raw <- read_excel("data/nipa_116a.xlsx")
gv_q_raw <- read_excel("data/nipa_116q.xlsx")
pop_raw <- read_excel("data/nipa_21q.xlsx")
gss_raw <- read_excel("data/gss_educ.xlsx")Confidence in Education
This document aims to explores the relationship between confidence in the US education system and actual government spending on education.
Loading Packages and Data
In this section, we pull in data downloaded from the Bureau of Economic Analysis (BEA) website and from the General Social Survey (GSS).
Cleaning data
First, lets clean the data on government spending over time in the US. Unfortunately, the BEA’s data on inflation adjusted spending on education only goes back to 2007. To work around this, we can take the ratio of nominal education spending to nominal total government spending and multiply real total government spending to get an estimate of real total government spending on education. However, total government spending on education is still misleading since it is increasing with things like population growth. To fix this, we can calculate education spending per capita by dividing by the U.S. population. We create two dataframes for the analysis: one with yearly data that we can use to compare to the yearly GSS data, and another that is quarterly for time series analysis.
### BEA NIPA data
## Population to get per capita spending
pop <- pop_raw %>%
select(-1) %>%
slice(-(1:4)) %>%
filter(row_number() %in% c(1,2,45)) %>%
t() %>% as.data.frame() %>%
slice(-1) %>%
rename(year = V1, quarter = V2, popq = V3) %>%
fill(year, .direction = "down") %>%
mutate(yearq = paste0(year,quarter),
year = as.numeric(year),
popq = as.numeric(popq)) %>%
select(year, yearq, popq)
pop_y <- pop %>%
group_by(year) %>%
summarise(popy = mean(popq, na.rm = TRUE))
## Getting ratio of education spending to total
govt_educ_y <- gved_y_raw %>%
select(-1) %>%
slice(-(1:4)) %>%
filter(row_number() %in% c(1,3,31)) %>%
t() %>% as.data.frame() %>%
slice(-1) %>%
rename(year = V1, govt = V2, educ = V3) %>%
mutate(educ = as.numeric(educ),
govt = as.numeric(govt),
year = as.numeric(year)) %>%
mutate(educ_govt = educ / govt)
## Getting annual government spending
govt_y <- gv_y_raw %>%
select(-1) %>%
slice(-(1:4)) %>%
filter(row_number() %in% c(1,24)) %>%
t() %>% as.data.frame() %>%
slice(-1) %>%
rename(year = V1, govty = V2) %>%
mutate(year = as.numeric(year),
govty = as.numeric(govty)) %>%
left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
filter(!is.na(educ_govt)) %>%
mutate(gved_y = as.numeric(govty) * as.numeric(educ_govt)) %>%
inner_join(pop_y, by = "year") %>%
mutate(gved_pc_y = (gved_y * 1e6) / (popy * 1e3)) %>%
select(year, gved_pc_y)
## Quarterly government spending on education (for time series forecast)
govt_q <- gv_q_raw %>%
select(-1) %>%
slice(-(1:4)) %>%
filter(row_number() %in% c(1,2,24)) %>%
t() %>% as.data.frame() %>%
slice(-1) %>%
rename(year = V1, quarter = V2, govtq = V3) %>%
fill(year, .direction = "down") %>%
mutate(yearq = paste0(year,quarter),
year = as.numeric(year),
quarter = as.numeric(substr(quarter,2,2)),
govtq = as.numeric(govtq)) %>%
select(c(year,quarter,yearq,govtq)) %>%
left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
filter(!is.na(educ_govt)) %>%
mutate(gved_q = govtq * educ_govt) %>%
inner_join(pop, by = c("yearq","year")) %>%
mutate(gved_pc_q = (gved_q * 1e6) / (popq * 1e3),
yearq = as.yearqtr(yearq)) %>%
select(c(year,quarter,yearq, gved_pc_q))Regarding the variables from the GSS, we use the following surveys:
- The coneduc survey asks whether the respondent has a great deal of confidence, only some confidence, or hardly any confidence in the institution of education in the US.
- The nateduc survey asks if the respondent thinks the US government is spending too much, too little, or the right amount on education.
We create two dataframes, one that has the average response per year and another that has the number of responses per response per year.
### GSS Education variables
## Taking an average of responses per year
gss_num <- gss_raw %>%
mutate(nateduc_num = case_when(
nateduc == "TOO LITTLE" ~ 1,
nateduc == "ABOUT RIGHT" ~ 2,
nateduc == "TOO MUCH" ~ 3,
TRUE ~ NA_real_
),
coneduc_num = case_when(
coneduc == "HARDLY ANY" ~ 1,
coneduc == "ONLY SOME" ~ 2,
coneduc == "A GREAT DEAL" ~ 3,
TRUE ~ NA_real_
),
year = as.numeric(year)) %>%
group_by(year) %>%
summarise(mean_nateduc = mean(nateduc_num, na.rm = TRUE),
mean_coneduc = mean(coneduc_num, na.rm = TRUE),
n_obs = n()) %>%
ungroup() %>%
mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE))
## Counting how many responses of each type per year
gss_count <- gss_raw %>%
mutate(nateduc = ifelse(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH"), nateduc, NA),
coneduc = ifelse(coneduc %in% c("A GREAT DEAL","HARDLY ANY","ONLY SOME"), coneduc, NA)) %>%
select(-id_) %>%
mutate(year = as.numeric(year)) %>%
pivot_longer(
cols = c(nateduc, coneduc),
names_to = "question",
values_to = "response"
) %>%
group_by(year, question, response) %>%
summarise(n = n(), .groups = "drop")Analysis
First, lets look at how government spending on education per capita has evolved over time.
## NIPA Data
ggplot(data = govt_y, aes(x=year,y=gved_pc_y)) +
geom_line() +
labs(title="Real Government Spending on Education Per Capita",
x = "Year",
y = "Chained (2017) Dollars") +
theme_minimal()This chart shows that the government has been spending more per capita on education over time. There were some dips in the 1980s and in the early 2010s, but overall it has been increasing. However, spending has stagnated a bit since the early 2000s.
Now that we’ve seen how spending has changed over time, lets look at confidence in education. We can make line charts for the proportion of each response the coneduc question.
## GSS Data
gss_prop <- gss_count %>%
filter(!is.na(response)) %>%
group_by(year, question) %>%
mutate(prop = n / sum(n))
# Confidence in Education
coneduc_prop <- gss_prop %>%
filter(question == "coneduc")
ggplot(coneduc_prop, aes(x = year, y = prop, color = response)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "How much confidence do you have in the education\nsystem in the US?",
x = "Year",
y = "Proportion of Respondents",
color = "Response"
) +
theme_bw() +
theme(axis.text.x = element_text(angle = 0, vjust = 0.5))It appears that confidence in the education system has decreased over time. The proportion of respondents who have a great deal of confidence has decreased from around 40% to less than 20% and the proportion who have hardly any confidence has increased from less than 10% to around 25%.
We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypotheses are:
\(H_0:\) The distribution of responses to how much confidence in education is the same prior to 1990 vs after 1990.
\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.
# Choose two periods for comparison (example: early vs recent)
gss_test1 <- gss_raw %>%
filter(coneduc %in% c("HARDLY ANY","ONLY SOME","A GREAT DEAL")) %>%
mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
group_by(period, coneduc) %>%
summarise(n = n(), .groups = "drop") %>%
pivot_wider(names_from = coneduc, values_from = n)
gss_test1# A tibble: 2 × 4
period `A GREAT DEAL` `HARDLY ANY` `ONLY SOME`
<chr> <int> <int> <int>
1 Early 7076 2366 11313
2 Late 7206 5488 17115
# Run the chi-square test
chisq.test(gss_test1[ , -1])
Pearson's Chi-squared test
data: gss_test1[, -1]
X-squared = 831.81, df = 2, p-value < 2.2e-16
Since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.
What about feelings about education spending? Let’s make a similar line chart over time to explore this.
# Spending on Education
nateduc_prop <- gss_prop %>%
filter(question == "nateduc")
ggplot(nateduc_prop, aes(x = year, y = prop, color = response)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Is US spending on education too much, too little, or about right?",
x = "Year",
y = "Proportion of Respondents",
color = "Response"
) +
theme_bw() +
theme(axis.text.x = element_text(angle = 0, vjust = 0.5))Since the beginning of the survey, most respondents believed that the government was spending too little on education, and this proportion has increased significantly to almost 80% of respondents.
We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypothese are:
\(H_0:\) The distribution of responses to how much are we spending on education is the same prior to 1990 vs after 1990.
\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.
# Choose two periods for comparison (example: early vs recent)
gss_test2 <- gss_raw %>%
filter(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH")) %>%
mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
group_by(period, nateduc) %>%
summarise(n = n(), .groups = "drop") %>%
pivot_wider(names_from = nateduc, values_from = n)
gss_test2# A tibble: 2 × 4
period `ABOUT RIGHT` `TOO LITTLE` `TOO MUCH`
<chr> <int> <int> <int>
1 Early 6183 10351 1440
2 Late 4989 17102 1432
# Run the chi-square test
chisq.test(gss_test2[ , -1])
Pearson's Chi-squared test
data: gss_test2[, -1]
X-squared = 1064.8, df = 2, p-value < 2.2e-16
Again, since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.
Finally, lets plot the mean of these two variables over time to create an even more clear visualization of the decline in sentiment.
gss_long <- gss_num %>%
pivot_longer(
cols = c(mean_nateduc, mean_coneduc),
names_to = "variable",
values_to = "mean_value"
)
ggplot(gss_long, aes(x = year, y = mean_value, color = variable)) +
geom_line(size = 1.1) +
labs(
title = "Mean Education Attitudes Over Time",
x = "Year",
y = "Mean Value",
color = "Variable"
) +
theme_minimal()Now that we’ve looked at the variables of interest separately, lets investigate how spending and confidence in education relate to each other. The following charts show how opinions about the amount of education spending relate to actual education spending.
nateduc_nipa <- gss_prop %>%
filter(question == "nateduc") %>%
left_join(govt_y, by = "year")
ggplot(nateduc_nipa, aes(x = gved_pc_y, y = prop)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~ response) +
labs(
title = "Confidence in amount of Education Spending vs Actual Spending",
x = "Real Education Spending Per Capita",
y = "Proportion of Respondents"
) +
scale_y_continuous(labels = scales::percent) +
theme_bw()`geom_smooth()` using formula = 'y ~ x'
The negative slope for the respondents who believe the government spends about the right amount on education and the positive slope for the those who believe the government spends too little on education indicate that even when the government prioritize education more, people still believe that the government isn’t spending enough. This may indicate that despite some efforts to spend more, the government still isn’t spending as much on education as they should.
Doing the same with confidence in the education system gives similar results.
coneduc_nipa <- gss_prop %>%
filter(question == "coneduc") %>%
left_join(govt_y, by = "year")
ggplot(coneduc_nipa, aes(x = gved_pc_y, y = prop)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~ response) +
labs(
title = "Confidence in the US Education System vs Spending",
x = "Real Education Spending Per Capita",
y = "Proportion of Respondents"
) +
scale_y_continuous(labels = scales::percent) +
theme_bw()`geom_smooth()` using formula = 'y ~ x'
Finally, we can do the same thing but using the mean of the surveys.
gss_num_govt <- gss_num %>%
inner_join(govt_y, by = "year") %>%
slice(-1)
# Mean_nateduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_nateduc)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "National Education Spending Attitudes vs Education Spending",
x = "Govt Spending per Capita",
y = "Mean National Education Attitude"
) +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
# Mean_coneduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_coneduc)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Confidence in Education vs Education Spending",
x = "Govt Spending per Capita",
y = "Mean Confidence in Education"
) +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
All these plots serve to show that even as government spending on education increases, people still increasingly have lost trust in the education system in the US. Let’s explore how these trends may evolve in the future with time series forecasts.
ARIMA
In this section, we analyze and attempt to forecast the time series for education spending and confidence in education. We focus only on confidence in education rather than feelings about speding on education since it more closely relates to our overall topic of trust in american institutions.
Government Spending on Education
gv_ts <- ts(govt_q$gved_pc_q,
start = c(min(govt_q$year),min(govt_q$quarter)),
frequency = 4) Lets start with a classical decomposition.
gv_ts_decomp <- decompose(gv_ts)
autoplot(gv_ts_decomp)The trend component of the chart shows a clear upward trend over time that matches the data. The seasonal component is also showing clear seasonality with spikes that appear to occur at the beginning of each year, however the effect of the seasonality seems small compared to the trend. Finally, the remainder shows fluctuations around zero that aren’t very uniform. It doesn’t seem likely that the remainder is showing cyclical patterns like recessions since it doesn’t seem to follow any particular pattern, but perhaps it reflects other outside political pressures on spending that are difficult to predict.
Let’s plot the lag and autocorrelation function plots to investigate the patterns more closely.
gglagplot(gv_ts, do.lines=FALSE, lags = 9) +
xlab("Lags") + ylab("Yt") +
ggtitle("Lag Plot for education spending") +
theme(axis.text.x=element_text(angle=45, hjust=1))The linear relationships in the lag plots suggest high serial correlation among consecutive observations.
ggAcf(gv_ts, 50) +
ggtitle("Autocorrelation Function for Real\nGovernment Spending on Education Per Capita") +
theme_minimal()This plot shows clear autocorrleation in the data, however the seasonality isn’t as present. Lets confirm stationarity with an augmented Dickey Fuller test.
tseries::adf.test(gv_ts)
Augmented Dickey-Fuller Test
data: gv_ts
Dickey-Fuller = -2.4123, Lag order = 6, p-value = 0.4023
alternative hypothesis: stationary
Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our suspicions from the ACF plot.
We can try to correct this with differencing.
diff_1 <- diff(gv_ts)
ggAcf(diff_1,50) +
ggtitle("ACF of First-Order Differenced Series") +
theme_minimal()It appears that one round of differencing may not have been sufficient since there is still a decreasing pattern every 4 lags. Let’s try one more time to be safe.
diff_2 <- diff(gv_ts, differences = 2)
ggAcf(diff_2,50) +
ggtitle("ACF of Second-Order Differenced Series") +
theme_minimal()ggPacf(diff_2,50) +
ggtitle("PACF of the Differenced Series") +
theme_minimal()This looks a little bit better, perhaps we will need to test for both 1 and 2 lags in the Arima. With this, we can proceed with the ARIMA.
First lets find the optimal parameters. Based on the differenced ACF and PACF plot, I believe the correct range for q is 0-8 and for p is 0-4.
# Define parameter ranges
p_range <- 0:4
d_range <- 1:2
q_range <- 0:9
# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)
# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)
# Initialize index for matrix row
i <- 1
# Loop through combinations of ARIMA model parameters
for (d in d_range) {
for (q in q_range) {
for (p in p_range) {
# Fit ARIMA model with specified (p,d,q)
model <- Arima(gv_ts, order = c(p, d, q), include.drift = TRUE)
# Store model parameters and AIC/BIC/AICc values in matrix
results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
# Increment row index
i <- i + 1
}
}
}
# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")
# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)
# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
kable_styling(full_width = FALSE, position = "center") %>%
row_spec(highlight_row, bold = TRUE, background = "#FFFF99") # Highlight row in yellow| p | d | q | AIC | BIC | AICc |
|---|---|---|---|---|---|
| 0 | 1 | 0 | 2491.868 | 2499.012 | 2491.914 |
| 1 | 1 | 0 | 2489.464 | 2500.181 | 2489.557 |
| 2 | 1 | 0 | 2491.169 | 2505.458 | 2491.324 |
| 3 | 1 | 0 | 2493.105 | 2510.965 | 2493.338 |
| 4 | 1 | 0 | 2413.508 | 2434.941 | 2413.836 |
| 0 | 1 | 1 | 2489.862 | 2500.578 | 2489.955 |
| 1 | 1 | 1 | 2477.350 | 2491.639 | 2477.506 |
| 2 | 1 | 1 | 2475.754 | 2493.614 | 2475.987 |
| 3 | 1 | 1 | 2469.481 | 2490.914 | 2469.809 |
| 4 | 1 | 1 | 2413.555 | 2438.560 | 2413.995 |
| 0 | 1 | 2 | 2491.114 | 2505.402 | 2491.269 |
| 1 | 1 | 2 | 2461.185 | 2479.046 | 2461.419 |
| 2 | 1 | 2 | 2471.032 | 2492.465 | 2471.361 |
| 3 | 1 | 2 | 2450.991 | 2475.996 | 2451.430 |
| 4 | 1 | 2 | 2415.378 | 2443.955 | 2415.945 |
| 0 | 1 | 3 | 2479.648 | 2497.509 | 2479.881 |
| 1 | 1 | 3 | 2453.874 | 2475.306 | 2454.202 |
| 2 | 1 | 3 | 2453.964 | 2478.969 | 2454.403 |
| 3 | 1 | 3 | 2448.673 | 2477.250 | 2449.240 |
| 4 | 1 | 3 | 2416.449 | 2448.599 | 2417.161 |
| 0 | 1 | 4 | 2426.689 | 2448.122 | 2427.017 |
| 1 | 1 | 4 | 2428.569 | 2453.574 | 2429.009 |
| 2 | 1 | 4 | 2427.455 | 2456.032 | 2428.022 |
| 3 | 1 | 4 | 2429.022 | 2461.172 | 2429.734 |
| 4 | 1 | 4 | 2418.418 | 2454.140 | 2419.291 |
| 0 | 1 | 5 | 2428.602 | 2453.607 | 2429.041 |
| 1 | 1 | 5 | 2424.873 | 2453.450 | 2425.440 |
| 2 | 1 | 5 | 2426.905 | 2459.054 | 2427.616 |
| 3 | 1 | 5 | 2427.978 | 2463.700 | 2428.851 |
| 4 | 1 | 5 | 2419.521 | 2458.815 | 2420.573 |
| 0 | 1 | 6 | 2429.563 | 2458.141 | 2430.130 |
| 1 | 1 | 6 | 2426.622 | 2458.771 | 2427.333 |
| 2 | 1 | 6 | 2423.968 | 2459.690 | 2424.841 |
| 3 | 1 | 6 | 2425.868 | 2465.161 | 2426.919 |
| 4 | 1 | 6 | 2421.496 | 2464.362 | 2422.744 |
| 0 | 1 | 7 | 2430.482 | 2462.632 | 2431.194 |
| 1 | 1 | 7 | 2428.499 | 2464.220 | 2429.372 |
| 2 | 1 | 7 | 2425.761 | 2465.055 | 2426.813 |
| 3 | 1 | 7 | 2427.382 | 2470.248 | 2428.630 |
| 4 | 1 | 7 | 2423.327 | 2469.765 | 2424.789 |
| 0 | 1 | 8 | 2426.343 | 2462.065 | 2427.216 |
| 1 | 1 | 8 | 2423.969 | 2463.263 | 2425.021 |
| 2 | 1 | 8 | 2425.870 | 2468.736 | 2427.118 |
| 3 | 1 | 8 | 2420.360 | 2466.798 | 2421.822 |
| 4 | 1 | 8 | 2423.777 | 2473.787 | 2425.471 |
| 0 | 1 | 9 | 2425.562 | 2464.856 | 2426.614 |
| 1 | 1 | 9 | 2425.747 | 2468.613 | 2426.995 |
| 2 | 1 | 9 | 2428.923 | 2475.361 | 2430.385 |
| 3 | 1 | 9 | 2422.604 | 2472.615 | 2424.298 |
| 4 | 1 | 9 | 2420.926 | 2474.508 | 2422.869 |
| 0 | 2 | 0 | 2694.714 | 2698.283 | 2694.730 |
| 1 | 2 | 0 | 2589.678 | 2596.815 | 2589.724 |
| 2 | 2 | 0 | 2556.299 | 2567.004 | 2556.392 |
| 3 | 2 | 0 | 2421.927 | 2436.201 | 2422.083 |
| 4 | 2 | 0 | 2418.246 | 2436.088 | 2418.481 |
| 0 | 2 | 1 | 2483.934 | 2491.071 | 2483.980 |
| 1 | 2 | 1 | 2471.136 | 2481.841 | 2471.229 |
| 2 | 2 | 1 | 2466.091 | 2480.365 | 2466.247 |
| 3 | 2 | 1 | 2416.039 | 2433.881 | 2416.274 |
| 4 | 2 | 1 | 2409.633 | 2431.043 | 2409.962 |
| 0 | 2 | 2 | 2461.556 | 2472.261 | 2461.649 |
| 1 | 2 | 2 | 2462.294 | 2476.568 | 2462.450 |
| 2 | 2 | 2 | 2458.668 | 2476.510 | 2458.902 |
| 3 | 2 | 2 | 2416.269 | 2437.680 | 2416.599 |
| 4 | 2 | 2 | 2409.464 | 2434.442 | 2409.905 |
| 0 | 2 | 3 | 2461.613 | 2475.886 | 2461.769 |
| 1 | 2 | 3 | 2463.498 | 2481.340 | 2463.732 |
| 2 | 2 | 3 | 2448.353 | 2469.763 | 2448.683 |
| 3 | 2 | 3 | 2418.269 | 2443.248 | 2418.710 |
| 4 | 2 | 3 | 2411.237 | 2439.784 | 2411.807 |
| 0 | 2 | 4 | 2460.940 | 2478.781 | 2461.174 |
| 1 | 2 | 4 | 2444.143 | 2465.554 | 2444.473 |
| 2 | 2 | 4 | 2446.544 | 2471.523 | 2446.985 |
| 3 | 2 | 4 | 2418.292 | 2446.839 | 2418.861 |
| 4 | 2 | 4 | 2412.406 | 2444.521 | 2413.121 |
| 0 | 2 | 5 | 2421.724 | 2443.134 | 2422.053 |
| 1 | 2 | 5 | 2422.362 | 2447.340 | 2422.803 |
| 2 | 2 | 5 | 2424.094 | 2452.640 | 2424.663 |
| 3 | 2 | 5 | 2420.282 | 2452.397 | 2420.996 |
| 4 | 2 | 5 | 2414.341 | 2450.025 | 2415.218 |
| 0 | 2 | 6 | 2422.448 | 2447.427 | 2422.889 |
| 1 | 2 | 6 | 2418.950 | 2447.497 | 2419.519 |
| 2 | 2 | 6 | 2420.950 | 2453.065 | 2421.664 |
| 3 | 2 | 6 | 2420.879 | 2456.562 | 2421.755 |
| 4 | 2 | 6 | 2415.589 | 2454.841 | 2416.645 |
| 0 | 2 | 7 | 2424.378 | 2452.924 | 2424.947 |
| 1 | 2 | 7 | 2420.950 | 2453.065 | 2421.664 |
| 2 | 2 | 7 | 2423.793 | 2459.477 | 2424.670 |
| 3 | 2 | 7 | 2420.419 | 2459.671 | 2421.475 |
| 4 | 2 | 7 | 2417.589 | 2460.409 | 2418.842 |
| 0 | 2 | 8 | 2426.250 | 2458.365 | 2426.964 |
| 1 | 2 | 8 | 2422.027 | 2457.710 | 2422.903 |
| 2 | 2 | 8 | 2422.764 | 2462.015 | 2423.820 |
| 3 | 2 | 8 | 2414.724 | 2457.545 | 2415.978 |
| 4 | 2 | 8 | 2419.284 | 2465.673 | 2420.752 |
| 0 | 2 | 9 | 2422.732 | 2458.416 | 2423.609 |
| 1 | 2 | 9 | 2420.357 | 2459.609 | 2421.413 |
| 2 | 2 | 9 | 2423.369 | 2466.189 | 2424.622 |
| 3 | 2 | 9 | 2415.815 | 2462.204 | 2417.283 |
| 4 | 2 | 9 | 2419.702 | 2469.659 | 2421.402 |
auto.arima(gv_ts)Series: gv_ts
ARIMA(1,2,1)(0,0,1)[4]
Coefficients:
ar1 ma1 sma1
-0.2286 -0.9327 0.4306
s.e. 0.0642 0.0314 0.0518
sigma^2 = 581.6: log likelihood = -1205.54
AIC=2419.09 AICc=2419.24 BIC=2433.36
The parameter search suggests an ARIMA(4,2,2) while the auto.arima() suggests ARIMA(1,2,1). Lets compare these to determine which one would be better.
model_output1 <- capture.output(sarima(gv_ts,4,2,2))# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1) # Locate where coefficient details start
end_line <- length(model_output1) # Last line of output
# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")Coefficients:
Estimate SE t.value p.value
ar1 0.0045 0.1008 0.0451 0.9641
ar2 0.0385 0.0536 0.7171 0.4740
ar3 0.0810 0.0530 1.5302 0.1272
ar4 0.5197 0.0540 9.6317 0.0000
ma1 -1.1826 0.1191 -9.9301 0.0000
ma2 0.1826 0.1182 1.5449 0.1236
sigma^2 estimated as 536.3868 on 256 degrees of freedom
AIC = 9.196427 AICc = 9.197684 BIC = 9.291764
model_output2 <- capture.output(sarima(gv_ts,1,2,1))# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output2) # Locate where coefficient details start
end_line <- length(model_output2) # Last line of output
# Print the relevant section automatically
cat(model_output2[start_line:end_line], sep = "\n")Coefficients:
Estimate SE t.value p.value
ar1 -0.2596 0.0652 -3.9825 1e-04
ma1 -0.8792 0.0360 -24.4325 0e+00
sigma^2 estimated as 708.7791 on 260 degrees of freedom
AIC = 9.431818 AICc = 9.431995 BIC = 9.472677
Both Residual Plots show nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.
Both Autocorrelation Functions (ACF) of the residuals reveal very few significant autocorrelations, however the plot for the ARIMA(1,2,1) does have a few more autocorrelations than the ARIMA(4,2,2).
Both Q-Q Plots indicate that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.
The Ljung-Box Test p-values for the ARIMA(4,2,2) remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals while in the ARIMA(1,2,1), all p-values remain below the 0.05 significance level implying some significant autocorrelations left in the residuals and concluding that model improvements may necessary.
Coefficient Significance: The model coefficients in the ARIMA(1,2,1) model are all significant while the coefficients for the first three AR terms and the last MA term are insignificant in the ARIMA(4,2,2) model.
fit1 <- Arima(gv_ts, order = c(4,2,2), include.drift = FALSE)
summary(fit1)Series: gv_ts
ARIMA(4,2,2)
Coefficients:
ar1 ar2 ar3 ar4 ma1 ma2
0.0045 0.0385 0.081 0.5197 -1.1826 0.1826
s.e. 0.1008 0.0536 0.053 0.0540 0.1191 0.1182
sigma^2 = 549: log likelihood = -1197.73
AIC=2409.46 AICc=2409.9 BIC=2434.44
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -1.015172 23.07227 17.61682 -0.03459217 0.8088008 0.3449994
ACF1
Training set -0.002555783
fit2 <- Arima(gv_ts, order = c(1,2,1), include.drift = FALSE)
summary(fit2)Series: gv_ts
ARIMA(1,2,1)
Coefficients:
ar1 ma1
-0.2596 -0.8792
s.e. 0.0652 0.0360
sigma^2 = 714.2: log likelihood = -1232.57
AIC=2471.14 AICc=2471.23 BIC=2481.84
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.2117617 26.52199 20.79503 0.01694402 0.9670959 0.4072402
ACF1
Training set -0.03802343
Looking at the error measurements, both the MAE and the RMSE are smaller in the ARIMA(4,2,2) model.
naive <- naive(gv_ts,h=5)
mean <- meanf(gv_ts,h=5)
drift <- rwf(gv_ts,drift=TRUE,h=5)
accuracy(naive) ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.835734 28.79736 21.69429 0.4615319 1.017407 0.4248508 -0.1290414
accuracy(mean) ME RMSE MAE MPE MAPE MASE
Training set -9.354218e-14 661.1524 594.0957 -11.0776 31.14832 11.63449
ACF1
Training set 0.9869592
accuracy(drift) ME RMSE MAE MPE MAPE MASE
Training set -1.085001e-13 27.40835 21.14088 0.04062096 0.9861303 0.4140131
ACF1
Training set -0.1290414
Furthermore, both models beat the benchmark models. Overall, it seems like the ARIMA(4,2,2) model is a better fit.
Now we can finally generate our forecast with the ARIMA(4,2,2) model.
# Generate the forecast
forecast_result1 <- forecast(fit1, h = 20)
# Plot the forecast
autoplot(forecast_result1) +
labs(title = "ARIMA(4,2,2) Forecast",
x = "Time",
y = "Real Education Spending per Capita") +
theme_minimal()This forecast suggests that real education spending per capita will continue to increase over time.
Educational Confidence
Now lets analyze confidence in the education system, since this is our primary variable of interest.
coneduc_ts_df <- gss_num %>%
select(year, mean_coneduc) %>%
arrange(year) %>%
complete(year = full_seq(year, 1)) %>%
mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE)) %>%
slice(-1)
coneduc_ts <- ts(coneduc_ts_df$mean_coneduc,
start = min(coneduc_ts_df$year),
end = max(coneduc_ts_df$year),
frequency = 1)Since this data is annual, we can’t use the decompose() function for a classical decomposition. Instead, we will extract the trend and remainder parts separately.
### Trend
trend <- stats::lowess(coneduc_ts)$y
remainder <- coneduc_ts - trend
# Plots
plot(coneduc_ts, type = "l",
main = "Data",
xlab = "Time", ylab = "Confidence in Education")plot(trend, type = "l",
main = "Trend Component",
xlab = "Time", ylab = "Trend")plot(remainder, type = "l",
main = "Remainder Component",
xlab = "Time", ylab = "Remainder")The trend component shows a clear downward trend in sentiment and the remainder shows some fluctuations around zero with some small spikes in the early years of the survey.
gglagplot(coneduc_ts, do.lines=FALSE, lags = 9) +
xlab("Lags") + ylab("Yt") +
ggtitle("Lag Plot for confidence in education") +
theme(axis.text.x=element_text(angle=45, hjust=1))The lag plots show a slight relationship between consecutive observations, suggesting there may not be very much autocorrelation.
ggAcf(coneduc_ts, 50) +
ggtitle("Autocorrelation Function for\nAverage Confidence in the Education System") +
theme_minimal()The ACF plot, however, shows some autocorrleation. Lets confirm stationarity with an augmented Dickey Fuller test.
tseries::adf.test(coneduc_ts)
Augmented Dickey-Fuller Test
data: coneduc_ts
Dickey-Fuller = -2.5281, Lag order = 3, p-value = 0.362
alternative hypothesis: stationary
Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our conclusions from the ACF plot.
Lets try to correct this with differencing.
diff_1 <- diff(coneduc_ts)
ggAcf(diff_1,50) +
ggtitle("ACF of First-Order Differenced Series") +
theme_minimal()ggPacf(diff_1,50) +
ggtitle("PACF of First-Order Differenced Series") +
theme_minimal()The ACF and PACF plots for the differenced series indicate that the data is now stationary and does not need additional differencing.
Now that we know that one difference is necessary, lets move on to a parameter search.
# Define parameter ranges
p_range <- 0:6
d_range <- 1
q_range <- 0:9
# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)
# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)
# Initialize index for matrix row
i <- 1
# Loop through combinations of ARIMA model parameters
for (d in d_range) {
for (q in q_range) {
for (p in p_range) {
# Fit ARIMA model with specified (p,d,q)
model <- Arima(coneduc_ts, order = c(p, d, q), include.drift = TRUE)
# Store model parameters and AIC/BIC/AICc values in matrix
results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
# Increment row index
i <- i + 1
}
}
}
# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")
# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)
# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
kable_styling(full_width = FALSE, position = "center") %>%
row_spec(highlight_row, bold = TRUE, background = "#FFFF99") # Highlight row in yellow| p | d | q | AIC | BIC | AICc |
|---|---|---|---|---|---|
| 0 | 1 | 0 | -146.1814 | -142.3177 | -145.9314 |
| 1 | 1 | 0 | -154.1756 | -148.3801 | -153.6650 |
| 2 | 1 | 0 | -160.2939 | -152.5666 | -159.4243 |
| 3 | 1 | 0 | -161.7735 | -152.1144 | -160.4402 |
| 4 | 1 | 0 | -160.6396 | -149.0486 | -158.7305 |
| 5 | 1 | 0 | -159.4146 | -145.8918 | -156.8099 |
| 6 | 1 | 0 | -161.3818 | -145.9272 | -157.9532 |
| 0 | 1 | 1 | -157.4889 | -151.6935 | -156.9783 |
| 1 | 1 | 1 | -155.6315 | -147.9042 | -154.7619 |
| 2 | 1 | 1 | -164.4317 | -154.7726 | -163.0984 |
| 3 | 1 | 1 | -162.4820 | -150.8911 | -160.5729 |
| 4 | 1 | 1 | -160.5157 | -146.9929 | -157.9110 |
| 5 | 1 | 1 | -159.1041 | -143.6495 | -155.6755 |
| 6 | 1 | 1 | -159.4157 | -142.0293 | -155.0255 |
| 0 | 1 | 2 | -155.9474 | -148.2201 | -155.0778 |
| 1 | 1 | 2 | -158.4396 | -148.7804 | -157.1062 |
| 2 | 1 | 2 | -162.4920 | -150.9010 | -160.5829 |
| 3 | 1 | 2 | -162.4265 | -148.9037 | -159.8218 |
| 4 | 1 | 2 | -160.9624 | -145.5078 | -157.5338 |
| 5 | 1 | 2 | -159.1660 | -141.7795 | -154.7757 |
| 6 | 1 | 2 | -159.7447 | -140.4264 | -154.2447 |
| 0 | 1 | 3 | -158.2586 | -148.5995 | -156.9253 |
| 1 | 1 | 3 | -159.5933 | -148.0023 | -157.6842 |
| 2 | 1 | 3 | -162.3082 | -148.7854 | -159.7035 |
| 3 | 1 | 3 | -161.0628 | -145.6082 | -157.6342 |
| 4 | 1 | 3 | -158.0646 | -140.6782 | -153.6744 |
| 5 | 1 | 3 | -156.1666 | -136.8484 | -150.6666 |
| 6 | 1 | 3 | -159.5780 | -138.3280 | -152.8088 |
| 0 | 1 | 4 | -162.1637 | -150.5727 | -160.2546 |
| 1 | 1 | 4 | -162.4235 | -148.9007 | -159.8188 |
| 2 | 1 | 4 | -160.5475 | -145.0929 | -157.1189 |
| 3 | 1 | 4 | -159.1997 | -141.8133 | -154.8095 |
| 4 | 1 | 4 | -158.1940 | -138.8757 | -152.6940 |
| 5 | 1 | 4 | -156.2203 | -134.9702 | -149.4510 |
| 6 | 1 | 4 | -160.4167 | -137.2348 | -152.2062 |
| 0 | 1 | 5 | -162.2605 | -148.7377 | -159.6559 |
| 1 | 1 | 5 | -161.2509 | -145.7963 | -157.8224 |
| 2 | 1 | 5 | -159.2674 | -141.8809 | -154.8771 |
| 3 | 1 | 5 | -159.8130 | -140.4947 | -154.3130 |
| 4 | 1 | 5 | -158.7977 | -137.5476 | -152.0284 |
| 5 | 1 | 5 | -157.0933 | -133.9114 | -148.8828 |
| 6 | 1 | 5 | -155.3629 | -130.2491 | -145.5250 |
| 0 | 1 | 6 | -160.3046 | -144.8500 | -156.8760 |
| 1 | 1 | 6 | -159.2570 | -141.8705 | -154.8667 |
| 2 | 1 | 6 | -157.4700 | -138.1517 | -151.9700 |
| 3 | 1 | 6 | -156.3756 | -135.1255 | -149.6063 |
| 4 | 1 | 6 | -157.4528 | -134.2709 | -149.2422 |
| 5 | 1 | 6 | -156.2534 | -131.1397 | -146.4156 |
| 6 | 1 | 6 | -157.7621 | -130.7165 | -146.0954 |
| 0 | 1 | 7 | -160.7412 | -143.3547 | -156.3509 |
| 1 | 1 | 7 | -158.7809 | -139.4626 | -153.2809 |
| 2 | 1 | 7 | -160.6976 | -139.4475 | -153.9284 |
| 3 | 1 | 7 | -158.6980 | -135.5161 | -150.4875 |
| 4 | 1 | 7 | -156.3615 | -131.2478 | -146.5237 |
| 5 | 1 | 7 | -154.9472 | -127.9016 | -143.2805 |
| 6 | 1 | 7 | -155.7651 | -126.7877 | -142.0508 |
| 0 | 1 | 8 | -158.8138 | -139.4956 | -153.3138 |
| 1 | 1 | 8 | -159.0816 | -137.8316 | -152.3124 |
| 2 | 1 | 8 | -158.6981 | -135.5162 | -150.4875 |
| 3 | 1 | 8 | -156.6980 | -131.5842 | -146.8601 |
| 4 | 1 | 8 | -157.0020 | -129.9564 | -145.3353 |
| 5 | 1 | 8 | -152.7827 | -123.8053 | -139.0684 |
| 6 | 1 | 8 | -154.9958 | -124.0866 | -138.9958 |
| 0 | 1 | 9 | -161.5029 | -140.2528 | -154.7337 |
| 1 | 1 | 9 | -159.6613 | -136.4794 | -151.4508 |
| 2 | 1 | 9 | -158.4304 | -133.3166 | -148.5925 |
| 3 | 1 | 9 | -152.4592 | -125.4136 | -140.7925 |
| 4 | 1 | 9 | -155.0022 | -126.0249 | -141.2880 |
| 5 | 1 | 9 | -153.0854 | -122.1762 | -137.0854 |
| 6 | 1 | 9 | -153.0750 | -120.2339 | -134.5295 |
auto.arima(coneduc_ts)Series: coneduc_ts
ARIMA(2,1,1) with drift
Coefficients:
ar1 ar2 ma1 drift
-1.0841 -0.6679 0.6965 -0.0074
s.e. 0.1485 0.1162 0.1436 0.0038
sigma^2 = 0.002027: log likelihood = 87.22
AIC=-164.43 AICc=-163.1 BIC=-154.77
According to the parameter search and the auto.arima() function, our best model is a simple ARIMA(2,1,1). Lets quickly look at the model diagnostics.
model_output1 <- capture.output(sarima(coneduc_ts,2,1,1))# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1) # Locate where coefficient details start
end_line <- length(model_output1) # Last line of output
# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")Coefficients:
Estimate SE t.value p.value
ar1 -1.0841 0.1485 -7.3013 0.0000
ar2 -0.6679 0.1162 -5.7458 0.0000
ma1 0.6965 0.1436 4.8502 0.0000
constant -0.0074 0.0038 -1.9609 0.0558
sigma^2 estimated as 0.001867907 on 47 degrees of freedom
AIC = -3.224151 AICc = -3.207101 BIC = -3.034756
The Residual Plot shows nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.
The Autocorrelation Function (ACF) of the residuals reveals no significant autocorrelations.
The Q-Q Plot indicates that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.
The Ljung-Box Test p-values remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals.
Coefficient Significance: The model coefficients in the ARIMA(2,1,1) model are all significant
Overall, this seems like a strong model. Lets compare it to some benchmark methods.
fit <- Arima(coneduc_ts, order = c(2,1,1), include.drift = FALSE)
summary(fit)Series: coneduc_ts
ARIMA(2,1,1)
Coefficients:
ar1 ar2 ma1
-1.0610 -0.6314 0.7172
s.e. 0.1504 0.1203 0.1375
sigma^2 = 0.002134: log likelihood = 85.43
AIC=-162.85 AICc=-161.98 BIC=-155.13
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.01067867 0.04438487 0.03073701 -0.5315693 1.437563 0.9150038
ACF1
Training set -0.01935951
naive <- naive(coneduc_ts,h=5)
mean <- meanf(coneduc_ts,h=5)
drift <- rwf(coneduc_ts,drift=TRUE,h=5)
accuracy(naive) ME RMSE MAE MPE MAPE MASE
Training set -0.006753008 0.05591304 0.03359222 -0.3509282 1.552349 1
ACF1
Training set -0.4017722
accuracy(mean) ME RMSE MAE MPE MAPE MASE
Training set -1.238466e-16 0.09197785 0.06769079 -0.1854871 3.184123 2.015073
ACF1
Training set 0.7513554
accuracy(drift) ME RMSE MAE MPE MAPE MASE
Training set 3.047601e-17 0.05550374 0.03343843 -0.03171947 1.541556 0.9954219
ACF1
Training set -0.4017722
Our ARIMA(2,1,1) beats all the benchmark methods according to the RMSE and the MAE.
Finally, we can plot our forecast.
# Generate the forecast
forecast_result <- forecast(fit, h = 5)
# Plot the forecast
autoplot(forecast_result) +
labs(title = "ARIMA(2,1,1) Forecast",
x = "Time",
y = "Average Confidence in Education") +
theme_minimal()The ARIMA predicts that confidence in education will remain low over time.
Conclusion
These analyses imply that that despite the government spending more on education per capita, people believe that education should be prioritized more. One possibility is that government spending on education is being misused and isn’t have as big of an effect on the quality of education as it should. Another possibility is that government spending in general hasn’t been increasing as much as it should so even if the portion of spending that is being used for education increases, the actual amount being spent isn’t meeting the demand.